
*******************************************************************************************
**project: 			Nighttime light as a proxy for development at the local level
**objective: 		Analysis of association between nighttime light and development indicators at the level of circular zones around DHS clusters and PRIO cells 
**author: 			AB
**date: 			May 2017
**last modified:	Feb 2018 (AB)
**use datasets:		DHSandLIGHT_PRIOcells.dta, DHSandLIGHT_BufferZones.dta
**output datasets: 	DHSandLIGHT_PRIOcells_overtime.dta
*******************************************************************************************

global path "D:\0. PUBLICATIONS\Journal articles\Empirical projects\PLOS 2018\Harvard Dataverse"
set matsize 1000
set more off, permanently
*set more on, permanently
cd "$path"

******************************************************************************************
**Table 1: Summary statistics
******************************************************************************************
use "$path/DHSandLIGHT_BufferZones.dta", clear
outreg2 using sumstats_bufferzones, sum(log) tex excel replace keep(wealth_clu NEwealth_clu NARprim_clu educationyears_clu IMR_3y_clu profbirth_clu MEAN_LIGHT MEAN_POPDENS electricity_clu urban)
**within-wave SD
bysort cw: egen sd_withincw = sd(MEAN_POPDENS)
bysort cw: gen cwindic=_n
sum sd_withincw if cwindic==1

use "$path/DHSandLIGHT_PRIOcells.dta", clear
outreg2 using sumstats_PRIOcells, sum(log) tex excel replace keep(wealth_cellmean NEwealth_cellmean NARprim_cellmean educationyears_cellmean IMR_3y_cellmean profbirth_cellmean MEAN_LIGHT MEAN_POPDENS electricity_cellmean urbanshare)
**within-wave SD
bysort cw: egen sd_withincw = sd(MEAN_POPDENS)
bysort cw: gen cwindic=_n
sum sd_withincw if cwindic==1

**List of countries and years
tab cw

/*
**Share of zero cluster means for each indicator
use "$path/DHSandLIGHT_BufferZones.dta", clear
local indic "NARprim educationyears IMR_3y profbirth housing wealth"
foreach ind of local indic {
	sum `ind'_clu
	local obs=`r(N)'
	count if `ind'_clu==0
	local zero=`r(N)'
	display `zero'/`obs'
	}
	
**Share of zero cellmean for each indicator
use "$path/DHSandLIGHT_PRIOcells.dta"
local indic "NARprim educationyears IMR_3y profbirth housing wealth"
foreach ind of local indic {
	sum `ind'_cellmean
	local obs=`r(N)' 
	count if `ind'_cellmean==0
	local zero=`r(N)'
	display `zero'/`obs'
	}		//35% zero-cells for IMR - better use winsorized indicators in main model
*/

******************************************************************************************
**Main results (Tables 2-7) 
******************************************************************************************

use "$path/DHSandLIGHT_BufferZones.dta", clear
gen IHS_LIGHT=ln(MEAN_LIGHT + ((MEAN_LIGHT^2+1)^0.5))
estimates clear
local ygroup "wealth NEwealth NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	qui: reghdfe `y'_clu ln_LIGHT01, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_1
	qui: reghdfe `y'_clu LIGHT_indic, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_2
	qui: reghdfe `y'_clu ln_LIGHT, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_3
	qui: reghdfe `y'_clu IHS_LIGHT, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_4
	qui: reghdfe `y'_clu ln_LIGHT01 ln_POPDENS, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_5
	qui: reghdfe `y'_clu ln_LIGHT01 ln_POPDENS electricity_clu urban, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_6
	estout Main_`y'_*, cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(r2 N) style(tex)
	}

use "$path/DHSandLIGHT_PRIOcells.dta", clear
gen IHS_LIGHT=ln(MEAN_LIGHT + ((MEAN_LIGHT^2+1)^0.5))
estimates clear
local ygroup "wealth NEwealth NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	qui: reghdfe `y'_cellmean ln_LIGHT01, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_1
	qui: reghdfe `y'_cellmean LIGHT_indic, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_2
	qui: reghdfe `y'_cellmean ln_LIGHT, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_3
	qui: reghdfe `y'_cellmean IHS_LIGHT, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_4
	qui: reghdfe `y'_cellmean ln_LIGHT01 ln_POPDENS, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_5
	qui: reghdfe `y'_cellmean ln_LIGHT01 ln_POPDENS electricity_cellmean urbanshare, a(cw_code) vce(cluster countryid year) 
	estimates store Main_`y'_6
	estout Main_`y'_*, cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(r2 N) style(tex)
	}
	
******************************************************************************************
**Robustness test: old satellites (Table S2)
******************************************************************************************
use "$path/DHSandLIGHT_BufferZones_oldsat.dta", clear
estimates clear
local ygroup "wealth NEwealth NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	qui: reghdfe `y'_clu ln_LIGHT01 ln_POPDENS, a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_1
	qui: reghdfe `y'_clu ln_LIGHT01 ln_POPDENS electricity_clu urban, a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_2
	}
estout R_*, cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(r2 N) style(tex)

use "$path/DHSandLIGHT_PRIOcells_oldsat.dta", clear
estimates clear
local ygroup "wealth NEwealth NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	qui: reghdfe `y'_cellmean ln_LIGHT01 ln_POPDENS, a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_1
	qui: reghdfe `y'_cellmean ln_LIGHT01 ln_POPDENS electricity_cellmean urbanshare, a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_2
	}
estout R_*, cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(r2 N) style(tex)

******************************************************************************************
**Robustness test: sum of lights (Table S3)
******************************************************************************************
use "$path/DHSandLIGHT_BufferZones.dta", clear
gen ln_LIGHTsum01=ln(COUNT_LIGHT*MEAN_LIGHT+.01)
cor ln_LIGHTsum01 ln_LIGHT01 
sum ln_LIGHTsum01 ln_LIGHT01 
estimates clear
local ygroup "wealth NEwealth NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	qui: reghdfe `y'_clu ln_LIGHTsum01 [aweight=overlap_inv], a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_1
	qui: reghdfe `y'_clu ln_LIGHTsum01 ln_POPDENS electricity_clu urban [aweight=overlap_inv], a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_2
	}
estout R_*, cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(r2 N) style(tex)

use "$path/DHSandLIGHT_PRIOcells.dta", clear
gen ln_LIGHTsum01=ln(COUNT_LIGHT*MEAN_LIGHT+.01)
cor ln_LIGHTsum01 ln_LIGHT01 
sum ln_LIGHTsum01 ln_LIGHT01 
estimates clear
local ygroup "wealth NEwealth NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	qui: reghdfe `y'_cellmean ln_LIGHTsum01 ln_POPDENS, a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_1
	qui: reghdfe `y'_cellmean ln_LIGHTsum01 ln_POPDENS electricity_cellmean urbanshare, a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_2
	}
estout R_*, cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(r2 N) style(tex)

******************************************************************************************
**Robustness test: cross-country (Table S4) 
******************************************************************************************
use "$path/DHSandLIGHT_BufferZones.dta", clear
estimates clear
local ygroup "NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	qui: reghdfe `y'_clu ln_LIGHT01 ln_POPDENS, a(year) vce(cluster countryid year) 
	estimates store R_`y'_1
	qui: reghdfe `y'_clu ln_LIGHT01 ln_POPDENS electricity_clu urban, a(year) vce(cluster countryid year) 
	estimates store R_`y'_2
	}
estout R_*, cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(r2 N) style(tex)

use "$path/DHSandLIGHT_PRIOcells.dta", clear
estimates clear
local ygroup "NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	qui: reghdfe `y'_cellmean ln_LIGHT01 ln_POPDENS, a(year) vce(cluster countryid year) 
	estimates store R_`y'_1
	qui: reghdfe `y'_cellmean ln_LIGHT01 ln_POPDENS electricity_cellmean urbanshare, a(year) vce(cluster countryid year) 
	estimates store R_`y'_2
	}
estout R_*, cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(r2 N) style(tex)

******************************************************************************************
**Robustness test: weighted by # observations  (Table S5)
******************************************************************************************
use "$path/DHSandLIGHT_BufferZones.dta", clear
estimates clear
local ygroup "wealth NEwealth NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	qui: reghdfe `y'_clu ln_LIGHT01 ln_POPDENS [aweight=count`y'_bas_clu], a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_1
	qui: reghdfe `y'_clu ln_LIGHT01 ln_POPDENS electricity_clu urban [aweight=count`y'_bas_clu], a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_2
	}
estout R_*, cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(r2 N) style(tex)

use "$path/DHSandLIGHT_PRIOcells.dta", clear
estimates clear
local ygroup "wealth NEwealth NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	qui: reghdfe `y'_cellmean ln_LIGHT01 ln_POPDENS [aweight=`y'_obsincell], a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_1
	qui: reghdfe `y'_cellmean ln_LIGHT01 ln_POPDENS electricity_cellmean urbanshare [aweight=`y'_obsincell], a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_2
	}
estout R_*, cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(r2 N) style(tex)

******************************************************************************************
**Robustness test: weighted by inverse # overlaps (Table S6, buffers only)
******************************************************************************************
use "$path/DHSandLIGHT_BufferZones.dta", clear
estimates clear
local ygroup "wealth NEwealth NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	qui: reghdfe `y'_clu ln_LIGHT01 ln_POPDENS [aweight=overlap_inv], a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_1
	qui: reghdfe `y'_clu ln_LIGHT01 ln_POPDENS electricity_clu urban [aweight=overlap_inv], a(cw_code) vce(cluster countryid year) 
	estimates store R_`y'_2
	}
estout R_*, cells(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.10 ** 0.05 *** 0.01) stats(r2 N) style(tex)

******************************************************************************************
**Figure 2: Distribution of development indicators within bins of light
******************************************************************************************
use "$path/DHSandLIGHT_BufferZones.dta", clear
xtile light_deciles = MEAN_LIGHT if MEAN_LIGHT>0, nquantiles(10)
tab light_deciles
bys light_deciles: sum MEAN_LIGHT
egen light_cat = cut(MEAN_LIGHT), at(0,0.017,0.637,2.216,5.09,9.17,16.01,26.056,37.321,49.85,60.4,63.1)
tab light_cat
local ygroup "NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	graph box `y'_clu, over(light_cat, relabel(1 "0" 2 "1" 3 "2" 4 "3" 5 "4" 6 "5" 7 "6" 8 "7" 9 "8" 10 "9" 11 "10") label(labsize(small))) ylabel(, labsize(small)) ytitle(, size(medium)) marker(1, msize(tiny)) saving(Fig_box_`y'_buffer, replace)
	graph export "$path/Fig_box_`y'_buffer.tif", replace
	graph export "$path/Fig_box_`y'_buffer.pdf", replace
	bysort light_cat: sum `y'_clu
	}

use "$path/DHSandLIGHT_PRIOcells.dta", clear
xtile light_deciles = MEAN_LIGHT if MEAN_LIGHT>0, nquantiles(10)
tab light_deciles
bys light_deciles: sum MEAN_LIGHT
egen light_cat = cut(MEAN_LIGHT), at(0, .0007,0.01158,0.0254,0.0445,0.0712,0.11192,0.187,0.3243,0.721,2.54,64)
tab light_cat
local ygroup "NARprim educationyears IMR_3y profbirth"
foreach y of local ygroup {
	graph box `y'_cellmean, alsize(100) over(light_cat, relabel(1 "0" 2 "1" 3 "2" 4 "3" 5 "4" 6 "5" 7 "6" 8 "7" 9 "8" 10 "9" 11 "10") label(labsize(small))) ylabel(, labsize(small)) ytitle(, size(medium)) marker(1, msize(tiny)) saving(Fig_box_`y'_cell, replace)
	graph export "$path/Fig_box_`y'_cell.tif", replace
	graph export "$path/Fig_box_`y'_cell.pdf", replace
	bysort light_cat: sum `y'_cellmean
	}	

graph combine "Fig_box_NARprim_buffer" "Fig_box_educationyears_buffer" "Fig_box_IMR_3y_buffer" "Fig_box_profbirth_buffer" "Fig_box_NARprim_cell" "Fig_box_educationyears_cell" "Fig_box_IMR_3y_cell" "Fig_box_profbirth_cell" , colfirst rows(4) cols(2) ysize(8.75) xsize(7)
graph export "$path/Fig_2.tif", width(11000) replace
graph export "$path/Fig_2.pdf", replace	

******************************************************************************************
**Figure 3: Effect heterogeneity wrt agriculture share in gdp
******************************************************************************************
use "$path/DHSandLIGHT_BufferZones.dta", clear
gen zvar=agrva
sum zvar
*for use in MVZ below: min 6.8, max 56.3
gen interaction=ln_LIGHT01*zvar
gen interaction_sq=interaction*zvar
reghdfe wealth_clu ln_LIGHT01 interaction interaction_sq ln_POPDENS, a(cw_code) vce(cluster countryid year) 
#delimit ;
matrix b=e(b);
matrix V=e(V);
scalar b1=b[1,1];
scalar b2=b[1,2];
scalar b3=b[1,3];
scalar varb1=V[1,1];
scalar varb2=V[2,2];
scalar varb3=V[3,3];
scalar covb1b2=V[1,2];
scalar covb1b3=V[1,3];
scalar covb2b3=V[2,3];
scalar list b1 b2 b3 varb1 varb2 varb3 covb1b2 covb1b3 covb2b3;
generate MVZ=(_n-1)/10;
replace  MVZ=. if _n>564;
replace  MVZ=. if _n<69;
gen conbx=b1+b2*MVZ+b3*(MVZ^2);
gen consx=sqrt(varb1+(varb2*(MVZ^2))+(varb3*(MVZ^4))+(2*covb1b2*MVZ)+(2*covb1b3*(MVZ^2))+(2*covb2b3*(MVZ^3)));
gen ax=1.96*consx;
gen upperx=conbx+ax;
gen lowerx=conbx-ax;
gen where=0;
egen tag_zvar = tag(zvar);
gen yline=0;
graph twoway hist zvar, width(2) percent color(gs14) yaxis(2)
	|| scatter where zvar if tag_zvar, plotr(m(zero)) ms(none) mlabcolor(gs5) mlabpos(6) legend(off)
	|| line conbx  MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1)
    || line upperx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| line lowerx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| ,
    xlabel(10 15 20 25 30 35 40 45 50 55, nogrid labsize(2))
    ylabel(0 0.1 0.2 0.3 0.4 0.5, axis(1) nogrid labsize(2))
    ylabel(0 5 10 15 20 25 30, axis(2) nogrid labsize(2))
	yscale(alt)
    yscale(alt axis(2))
    legend(off)
    xtitle("Agriculture share in GDP", size(3.5))
    ytitle("Marginal effect", size(3))
	ytitle("Observations in %", axis(2) size(3))
    xsca(titlegap(2))
    ysca(titlegap(2))
	saving(Fig_agr_buffer, replace);
graph export "$path/Fig_agr_buffer.tif", replace;
graph export "$path/Fig_agr_buffer.pdf", replace;

#delimit cr
		
use "$path/DHSandLIGHT_PRIOcells.dta", clear
gen zvar=agrva
sum zvar
*for use in MVZ below: min 6.8, max 56.3
gen interaction=ln_LIGHT01*zvar
gen interaction_sq=interaction*zvar
reghdfe wealth_cellmean ln_LIGHT01 interaction interaction_sq ln_POPDENS, a(cw_code) vce(cluster countryid year) 
#delimit ;
matrix b=e(b);
matrix V=e(V);
scalar b1=b[1,1];
scalar b2=b[1,2];
scalar b3=b[1,3];
scalar varb1=V[1,1];
scalar varb2=V[2,2];
scalar varb3=V[3,3];
scalar covb1b2=V[1,2];
scalar covb1b3=V[1,3];
scalar covb2b3=V[2,3];
scalar list b1 b2 b3 varb1 varb2 varb3 covb1b2 covb1b3 covb2b3;
generate MVZ=(_n-1)/10;
replace  MVZ=. if _n>564;
replace  MVZ=. if _n<69;
gen conbx=b1+b2*MVZ+b3*(MVZ^2);
gen consx=sqrt(varb1+(varb2*(MVZ^2))+(varb3*(MVZ^4))+(2*covb1b2*MVZ)+(2*covb1b3*(MVZ^2))+(2*covb2b3*(MVZ^3)));
gen ax=1.96*consx;
gen upperx=conbx+ax;
gen lowerx=conbx-ax;
gen where=0;
egen tag_zvar = tag(zvar);
gen yline=0;
graph twoway hist zvar, width(2) percent color(gs14) yaxis(2)
	|| scatter where zvar if tag_zvar, plotr(m(zero)) ms(none) mlabcolor(gs5) mlabpos(6) legend(off)
	|| line conbx  MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1)
    || line upperx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| line lowerx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| ,
    xlabel(10 15 20 25 30 35 40 45 50 55, nogrid labsize(2))
    ylabel(0 0.1 0.2 0.3 0.4 0.5, axis(1) nogrid labsize(2))
    ylabel(0 5 10 15 20 25 30, axis(2) nogrid labsize(2))
	yscale(alt)
    yscale(alt axis(2))
    legend(off)
    xtitle("Agriculture share in GDP", size(3.5))
    ytitle("Marginal effect", size(3))
	ytitle("Observations in %" , axis(2) size(3))
    xsca(titlegap(2))
    ysca(titlegap(2))
	saving(Fig_agr_cell, replace);
graph export "$path/Fig_agr_cell.tif", replace;
graph export "$path/Fig_agr_cell.pdf", replace;
#delimit cr

******************************************************************************************
**Figure 3: Effect heterogeneity wrt polity
******************************************************************************************	
use "$path/DHSandLIGHT_BufferZones.dta", clear
*pol is based on freedom house, which has transformed the original polity2 score.
*We need to re-engineer the original polity2 score
gen polity2=2*pol-10
gen zvar=polity2
sum zvar
gen interaction=ln_LIGHT01*zvar
gen interaction_sq=interaction*zvar
reghdfe wealth_clu ln_LIGHT01 interaction interaction_sq ln_POPDENS, a(cw_code) vce(cluster countryid year) 
#delimit ;
matrix b=e(b);
matrix V=e(V);
scalar b1=b[1,1];
scalar b2=b[1,2];
scalar b3=b[1,3];
scalar varb1=V[1,1];
scalar varb2=V[2,2];
scalar varb3=V[3,3];
scalar covb1b2=V[1,2];
scalar covb1b3=V[1,3];
scalar covb2b3=V[2,3];
scalar list b1 b2 b3 varb1 varb2 varb3 covb1b2 covb1b3 covb2b3;
generate MVZ=(_n-800)/100;
replace  MVZ=. if _n>1720;
replace  MVZ=. if _n<-800;
gen conbx=b1+b2*MVZ+b3*(MVZ^2);
gen consx=sqrt(varb1+(varb2*(MVZ^2))+(varb3*(MVZ^4))+(2*covb1b2*MVZ)+(2*covb1b3*(MVZ^2))+(2*covb2b3*(MVZ^3)));
gen ax=1.96*consx;
gen upperx=conbx+ax;
gen lowerx=conbx-ax;
gen where=0;
egen tag_zvar = tag(zvar);
gen yline=0;
graph twoway hist zvar, width(1) percent color(gs14) yaxis(2)
	|| scatter where zvar if tag_zvar, plotr(m(zero)) ms(none) mlabcolor(gs5) mlabpos(6) legend(off)
	|| line conbx  MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1)
    || line upperx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| line lowerx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| ,
    xlabel(-8 -6 -4 -2 0 2 4 6 8, nogrid labsize(2))
    ylabel(0 0.1 0.2 0.3 0.4 0.5, axis(1) nogrid labsize(2))
    ylabel(0 5 10 15 20 25 30, axis(2) nogrid labsize(2))
	yscale(alt)
    yscale(alt axis(2))
    legend(off)
    xtitle("Polity2 score", size(3.5))
    ytitle("Marginal effect", size(3))
	ytitle("Observations in %" , axis(2) size(3))
    xsca(titlegap(2))
    ysca(titlegap(2))
	saving(Fig_pol_buffer, replace);
graph export "$path/Fig_pol_buffer.tif", replace;
graph export "$path/Fig_pol_buffer.pdf", replace;
#delimit cr
	
use "$path/DHSandLIGHT_PRIOcells.dta", clear
gen polity2=2*pol-10
gen zvar=polity2
sum zvar
gen interaction=ln_LIGHT01*zvar
gen interaction_sq=interaction*zvar
reghdfe wealth_cellmean ln_LIGHT01 interaction interaction_sq ln_POPDENS, a(cw_code) vce(cluster countryid year) 
#delimit ;
matrix b=e(b);
matrix V=e(V);
scalar b1=b[1,1];
scalar b2=b[1,2];
scalar b3=b[1,3];
scalar varb1=V[1,1];
scalar varb2=V[2,2];
scalar varb3=V[3,3];
scalar covb1b2=V[1,2];
scalar covb1b3=V[1,3];
scalar covb2b3=V[2,3];
scalar list b1 b2 b3 varb1 varb2 varb3 covb1b2 covb1b3 covb2b3;
generate MVZ=(_n-800)/100;
replace  MVZ=. if _n>1720;
replace  MVZ=. if _n<-800;
gen conbx=b1+b2*MVZ+b3*(MVZ^2);
gen consx=sqrt(varb1+(varb2*(MVZ^2))+(varb3*(MVZ^4))+(2*covb1b2*MVZ)+(2*covb1b3*(MVZ^2))+(2*covb2b3*(MVZ^3)));
gen ax=1.96*consx;
gen upperx=conbx+ax;
gen lowerx=conbx-ax;
gen where=0;
egen tag_zvar = tag(zvar);
gen yline=0;
graph twoway hist zvar, width(1) percent color(gs14) yaxis(2)
	|| scatter where zvar if tag_zvar, plotr(m(zero)) ms(none) mlabcolor(gs5) mlabpos(6) legend(off)
	|| line conbx  MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1)
    || line upperx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| line lowerx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| ,
    xlabel(-8 -6 -4 -2 0 2 4 6 8, nogrid labsize(2))
    ylabel(0 0.1 0.2 0.3 0.4 0.5, axis(1) nogrid labsize(2))
    ylabel(0 5 10 15 20 25 30, axis(2) nogrid labsize(2))
	yscale(alt)
    yscale(alt axis(2))
    legend(off)
    xtitle("Polity2 score", size(3.5))
    ytitle("Marginal effect", size(3))
	ytitle("Observations in %" , axis(2) size(3))
    xsca(titlegap(2))
    ysca(titlegap(2))
	saving(Fig_pol_cell, replace);
graph export "$path/Fig_pol_cell.tif", replace;
graph export "$path/Fig_pol_cell.pdf", replace;
#delimit cr
		
******************************************************************************************
**Figure 3: Effect heterogeneity wrt log gdp
******************************************************************************************
use "$path/DHSandLIGHT_BufferZones.dta", clear
gen zvar=ln(gdp)
sum zvar
*for use in MVZ below: min 4.92, max 8.42
gen interaction=ln_LIGHT01*zvar
gen interaction_sq=interaction*zvar
reghdfe wealth_clu ln_LIGHT01 interaction interaction_sq ln_POPDENS, a(cw_code) vce(cluster countryid year) 
#delimit ;
matrix b=e(b);
matrix V=e(V);
scalar b1=b[1,1];
scalar b2=b[1,2];
scalar b3=b[1,3];
scalar varb1=V[1,1];
scalar varb2=V[2,2];
scalar varb3=V[3,3];
scalar covb1b2=V[1,2];
scalar covb1b3=V[1,3];
scalar covb2b3=V[2,3];
scalar list b1 b2 b3 varb1 varb2 varb3 covb1b2 covb1b3 covb2b3;
generate MVZ=(_n-1)/100;
replace  MVZ=. if _n>843;
replace  MVZ=. if _n<493;
gen conbx=b1+b2*MVZ+b3*(MVZ^2);
gen consx=sqrt(varb1+(varb2*(MVZ^2))+(varb3*(MVZ^4))+(2*covb1b2*MVZ)+(2*covb1b3*(MVZ^2))+(2*covb2b3*(MVZ^3)));
gen ax=1.96*consx;
gen upperx=conbx+ax;
gen lowerx=conbx-ax;
gen where=0;
egen tag_zvar = tag(zvar);
gen yline=0;
graph twoway hist zvar, width(0.1) percent color(gs14) yaxis(2)
	|| scatter where zvar if tag_zvar, plotr(m(zero)) ms(none) mlabcolor(gs5) mlabpos(6) legend(off)
	|| line conbx  MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1)
    || line upperx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| line lowerx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| ,
    xlabel(5.0 5.5 6.0 6.5 7.0 7.5 8.0, nogrid labsize(2))
    ylabel(0 0.1 0.2 0.3 0.4 0.5, axis(1) nogrid labsize(2))
    ylabel(0 5 10 15 20 25 30, axis(2) nogrid labsize(2))
	yscale(alt)
    yscale(alt axis(2))
    legend(off)
    xtitle("ln(GDP per capita)", size(3.5))
    ytitle("Marginal effect", size(3))
	ytitle("Observations in %" , axis(2) size(3))
    xsca(titlegap(2))
    ysca(titlegap(2))
	saving(Fig_gdp_buffer, replace);
graph export "$path/Fig_gdp_buffer.tif", replace;
graph export "$path/Fig_gdp_buffer.pdf", replace;
#delimit cr

use "$path/DHSandLIGHT_PRIOcells.dta", clear
gen zvar=ln(gdp)
sum zvar
*for use in MVZ below: min 4.917, max 8.412
gen interaction=ln_LIGHT01*zvar
gen interaction_sq=interaction*zvar
reghdfe wealth_cellmean ln_LIGHT01 interaction interaction_sq ln_POPDENS, a(cw_code) vce(cluster countryid year) 
#delimit ;
matrix b=e(b);
matrix V=e(V);
scalar b1=b[1,1];
scalar b2=b[1,2];
scalar b3=b[1,3];
scalar varb1=V[1,1];
scalar varb2=V[2,2];
scalar varb3=V[3,3];
scalar covb1b2=V[1,2];
scalar covb1b3=V[1,3];
scalar covb2b3=V[2,3];
scalar list b1 b2 b3 varb1 varb2 varb3 covb1b2 covb1b3 covb2b3;
generate MVZ=(_n-1)/100;
replace  MVZ=. if _n>843;
replace  MVZ=. if _n<493;
gen conbx=b1+b2*MVZ+b3*(MVZ^2);
gen consx=sqrt(varb1+(varb2*(MVZ^2))+(varb3*(MVZ^4))+(2*covb1b2*MVZ)+(2*covb1b3*(MVZ^2))+(2*covb2b3*(MVZ^3)));
gen ax=1.96*consx;
gen upperx=conbx+ax;
gen lowerx=conbx-ax;
gen where=0;
egen tag_zvar = tag(zvar);
gen yline=0;
graph twoway hist zvar, width(0.1) percent color(gs14) yaxis(2)
	|| scatter where zvar if tag_zvar, plotr(m(zero)) ms(none) mlabcolor(gs5) mlabpos(6) legend(off)
	|| line conbx  MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1)
    || line upperx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| line lowerx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| ,
    xlabel(5.0 5.5 6.0 6.5 7.0 7.5 8.0, nogrid labsize(2))
    ylabel(0 0.1 0.2 0.3 0.4 0.5, axis(1) nogrid labsize(2))
    ylabel(0 5 10 15 20 25 30, axis(2) nogrid labsize(2))
	yscale(alt)
    yscale(alt axis(2))
    legend(off)
    xtitle("ln(GDP per capita)", size(3.5))
    ytitle("Marginal effect", size(3))
	ytitle("Observations in %" , axis(2) size(3))
    xsca(titlegap(2))
    ysca(titlegap(2))
	saving(Fig_gdp_cell, replace);
graph export "$path/Fig_gdp_cell.tif", replace;
graph export "$path/Fig_gdp_cell.pdf", replace;
#delimit cr
	
******************************************************************************************
**Figure 3: Effect heterogeneity wrt year
******************************************************************************************	
use "$path/DHSandLIGHT_BufferZones.dta", clear
gen time=year-1992
gen zvar=time
gen interaction=ln_LIGHT01*zvar
gen interaction_sq=interaction*zvar
reghdfe wealth_clu ln_LIGHT01 interaction interaction_sq ln_POPDENS, a(cw_code) vce(cluster countryid year) 
#delimit ;
matrix b=e(b);
matrix V=e(V);
scalar b1=b[1,1];
scalar b2=b[1,2];
scalar b3=b[1,3];
scalar varb1=V[1,1];
scalar varb2=V[2,2];
scalar varb3=V[3,3];
scalar covb1b2=V[1,2];
scalar covb1b3=V[1,3];
scalar covb2b3=V[2,3];
scalar list b1 b2 b3 varb1 varb2 varb3 covb1b2 covb1b3 covb2b3;
generate MVZ=(_n-1)/100;
replace  MVZ=. if _n>2101;
replace  MVZ=. if _n<0;
gen conbx=b1+b2*MVZ+b3*(MVZ^2);
gen consx=sqrt(varb1+(varb2*(MVZ^2))+(varb3*(MVZ^4))+(2*covb1b2*MVZ)+(2*covb1b3*(MVZ^2))+(2*covb2b3*(MVZ^3)));
gen ax=1.96*consx;
gen upperx=conbx+ax;
gen lowerx=conbx-ax;
gen where=0;
egen tag_zvar = tag(zvar);
gen yline=0;
graph twoway hist zvar, width(1) percent color(gs14) yaxis(2)
	|| scatter where zvar if tag_zvar, plotr(m(zero)) ms(none) mlabcolor(gs5) mlabpos(6) legend(off)
	|| line conbx  MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1)
    || line upperx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| line lowerx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| ,
    xlabel(0 2 4 6 8 10 12 14 16 18 20, nogrid labsize(2))
	ylabel(0 0.1 0.2 0.3 0.4 0.5, axis(1) nogrid labsize(2))
    ylabel(0 5 10 15 20 25 30, axis(2) nogrid labsize(2))
 	yscale(alt)
    yscale(alt axis(2))
	legend(off)
    xtitle("Years since 1992", size(3.5))
    ytitle("Marginal effect", size(3))
	ytitle("Observations in %" , axis(2) size(3))
    xsca(titlegap(2))
    ysca(titlegap(2))
	saving(Fig_years_buffer, replace);
graph export "$path/Fig_years_buffer.tif", replace;
graph export "$path/Fig_years_buffer.pdf", replace;
#delimit cr
	
use "$path/DHSandLIGHT_PRIOcells.dta", clear
gen time=year-1992
gen zvar=time
gen interaction=ln_LIGHT01*zvar
gen interaction_sq=interaction*zvar
reghdfe wealth_cellmean ln_LIGHT01 interaction interaction_sq ln_POPDENS, a(cw_code) vce(cluster countryid year) 
#delimit ;
matrix b=e(b);
matrix V=e(V);
scalar b1=b[1,1];
scalar b2=b[1,2];
scalar b3=b[1,3];
scalar varb1=V[1,1];
scalar varb2=V[2,2];
scalar varb3=V[3,3];
scalar covb1b2=V[1,2];
scalar covb1b3=V[1,3];
scalar covb2b3=V[2,3];
scalar list b1 b2 b3 varb1 varb2 varb3 covb1b2 covb1b3 covb2b3;
generate MVZ=(_n-1)/100;
replace  MVZ=. if _n>2101;
replace  MVZ=. if _n<0;
gen conbx=b1+b2*MVZ+b3*(MVZ^2);
gen consx=sqrt(varb1+(varb2*(MVZ^2))+(varb3*(MVZ^4))+(2*covb1b2*MVZ)+(2*covb1b3*(MVZ^2))+(2*covb2b3*(MVZ^3)));
gen ax=1.96*consx;
gen upperx=conbx+ax;
gen lowerx=conbx-ax;
gen where=0;
egen tag_zvar = tag(zvar);
gen yline=0;
graph twoway hist zvar, width(1) percent color(gs14) yaxis(2)
	|| scatter where zvar if tag_zvar, plotr(m(zero)) ms(none) mlabcolor(gs5) mlabpos(6) legend(off)
	|| line conbx  MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1)
    || line upperx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| line lowerx MVZ, clpattern(dash) clwidth(thin) clcolor(black)
	|| ,
    xlabel(0 2 4 6 8 10 12 14 16 18 20, nogrid labsize(2))
    ylabel(0 0.1 0.2 0.3 0.4 0.5, axis(1) nogrid labsize(2))
    ylabel(0 5 10 15 20 25 30, axis(2) nogrid labsize(2))
 	yscale(alt)
    yscale(alt axis(2))
	legend(off)
    xtitle("Years since 1992", size(3.5))
    ytitle("Marginal effect", size(3))
	ytitle("Observations in %" , axis(2) size(3))
    xsca(titlegap(2))
    ysca(titlegap(2))
	saving(Fig_years_cell, replace);
graph export "$path/Fig_years_cell.tif", replace;
graph export "$path/Fig_years_cell.pdf", replace;
#delimit cr	

graph combine "Fig_gdp_buffer" "Fig_gdp_cell" "Fig_agr_buffer" "Fig_agr_cell" "Fig_pol_buffer" "Fig_pol_cell" "Fig_years_buffer" "Fig_years_cell", rows(1) cols(2) ysize(8.75) xsize(7)
graph export "$path/Fig_3.tif", width(11000) replace
graph export "$path/Fig_3.pdf", replace



******************************************************************************************
**Figure SM: Dark clusters and grid cells
******************************************************************************************

**Dark clusters over time
use "$path/DHSandLIGHT_BufferZones.dta", clear
replace LIGHT_indic=1-LIGHT_indic
sum LIGHT_indic
**35 surveys up to 2005, 36 from 2006
sum LIGHT_indic if year<2006
sum LIGHT_indic if year>2005
sum LIGHT_indic if year<2006 & cc!="eg"
sum LIGHT_indic if year>2005 & cc!="eg"

collapse LIGHT_indic, by(cc countryid year)
sum LIGHT_indic
#delimit ;
graph twoway scatter (LIGHT_indic year), msymbol(i) mlabposition(0) mlabel(cc)  
	|| lfit (LIGHT_indic year), clpattern(solid) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="eg", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="ml", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="bf", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="bj", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="ci", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="et", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="gh", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="gn", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="mw", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="ng", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="nm", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="sn", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="ug", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="zw", clpattern(dash) clwidth(thin) clcolor(black)
	|| ,
	xlabel(1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 2014, nogrid labsize(3))
    ylabel(0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9, nogrid labsize(3))
    legend(off)
    xtitle("Year", size(3.5))
    ytitle("Share of spatial units with light=0", size(3.5))
	scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white));
graph export "$path/Fig_dark_cluster.pdf", replace;
#delimit cr	


**Dark cells over time
use "$path/DHSandLIGHT_PRIOcells.dta", clear
replace LIGHT_indic=1-LIGHT_indic
sum LIGHT_indic
**35 surveys up to 2005, 36 from 2006
sum LIGHT_indic if year<2006
sum LIGHT_indic if year>2005
sum LIGHT_indic if year<2006 & cc!="eg"
sum LIGHT_indic if year>2005 & cc!="eg"

collapse LIGHT_indic, by(cc countryid year)
sum LIGHT_indic
#delimit ;
graph twoway scatter (LIGHT_indic year), msymbol(i) mlabposition(0) mlabel(cc)  
	|| lfit (LIGHT_indic year), clpattern(solid) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="eg", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="ml", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="bf", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="bj", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="ci", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="et", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="gh", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="gn", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="mw", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="ng", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="nm", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="sn", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="ug", clpattern(dash) clwidth(thin) clcolor(black)
	|| lfit (LIGHT_indic year) if cc=="zw", clpattern(dash) clwidth(thin) clcolor(black)
	|| ,
	xlabel(1992 1994 1996 1998 2000 2002 2004 2006 2008 2010 2012 2014, nogrid labsize(3))
    ylabel(0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9, nogrid labsize(3))
    legend(off)
    xtitle("Year", size(3.5))
    ytitle("Share of spatial units with light=0", size(3.5))
	scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white));
graph export "$path/Fig_dark_cell.pdf", replace;
#delimit cr
